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A Eulerian-Lagrangian model to simulate 
two-phase/particulate flows 

By S. V. Apte, K. Maheshf, & T. Lundgrenf 


1. Motivation and Objectives 

Figure 1 shows a snapshot of liquid fuel spray coming out of an injector nozzle in 
a realistic gas-turbine combustor. Here the spray atomization was simulated using a 
stochastic secondary breakup model (Apte et al. 2003a) with point-particle approxima- 
tion for the droplets. Very close to the injector, it is observed that the spray density is 
large and the droplets cannot be treated as point-particles. The volume displaced by the 
liquid in this region is significant and can alter the gas-phase flow and spray evolution. 
In order to address this issue, one can compute the dense spray regime by an Eulerian- 
Eulerian technique using advanced interface tracking/level-set methods (Sussman et al. 
1994; Tryggvason et al. 2001; Herrmann 2003). This, however, is computationally in- 
tensive and may not be viable in realistic complex configurations. We therefore plan to 
develop a methodology based on Eulerian-Lagrangian technique which will allow us to 
capture the essential features of primary atomization using models to capture interac- 
tions between the fluid and droplets and which can be directly applied to the standard 
atomization models used in practice. The numerical scheme for unstructured grids de- 
veloped by Mahesh et al. (2003) for incompressible flows is modified to take into account 
the droplet volume fraction. The numerical framework is directly applicable to realistic 
combustor geometries. 

Our main objectives in this work are: 

• Develop a numerical formulation based on Eulerian-Lagrangian techniques with 
models for interaction terms between the fluid and particles to capture the Kelvin- 
Helmholtz type instabilities observed during primary atomization. 

• Validate this technique for various two-phase and particulate flows. 

• Assess its applicability to capture primary atomization of liquid jets in conjuction 
with secondary atomization models. 


2. Mathematical Formulation 

Recent direct numerical simulations of large number of solid particles interacting 
through a fluid medium by Joseph and collaborators (Choi & Joseph 2001) show that a 
layer of heavy particles with fluid streaming above it can develop Kelvin-Helmholtz (K- 
H) instability waves whereas a layer of particles above a lighter fluid develops Rayleigh- 
Taylor instability. This suggests that the primary breakup of a liquid jet into a spray can 
be simulated by replacing the jet by a closely packed collection of droplets with some 
assumed size distribution. The K-H instability at the boundary between droplets and 
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Figure 1. Snapshot spray from a gas-turbine fuel-injector. 


fluid would initiate dispersal of droplets into a spray. Further breakup of these dispersed 
drops can be obtained by advanced secondary breakup models. 

The formulation described below is a modification of the equations for spray compu- 
tation developed by Dukowicz (1980) which consists of Eulerian fluid and Lagrangian 
particle calculations, and accounts for the displacement of the fluid by the particles as 
well as the momentum interchange between them. The modifications presented here are 
mainly in the details of modeling the interaction terms. 

2.1. Gas-Phase Equations 

The fluid mass for unit volume satisfies a continuity equation, 

(P/©/) + V ’ (P/@/ u /) = 0 (2-1) 

where pf, Qf , and u / are fluid density, volume fraction, and velocity, respectively. This 
indicates that the average velocity field of the fluid phase does not satisfy the divergence- 
free condition even if we consider an incompressible suspending fluid. The fluid momen- 
tum equation is given as 

^ (p/©/u/) + V • (pfQfUfUf) = - V (©/p) + V • (P/Dc) + F (2.2) 

where p is the average dynamic pressure in the fluid phase, pf is the viscosity of the 
fluid phase, and D c = V u c + V u ^ is the average deformation-rate of the fluid-particle 
composite, u c is the composite velocity of the mixture, and F is the force per unit 
volume the particles exert on gas. These equations are derived in detail for constant 
density flows by Joseph & Lundgren (1990). For particulate flows and dilute suspensions 
at low Reynolds numbers, the fluid viscosity should be replaced by an effective viscosity 
p* by using Thomas (1965) correlation, 

p* =p f ( 1 + 2.50/ + 10.050/ + 0.00273e 16 ' 6O/ ) (2.3) 

2.2. Particle-Phase Equations 

The evolution of particle-phase is governed by a Liouville equation for the particle dis- 
tribution function <f>(x p , u p ,p p ,V p ,t) 
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— + Vx-($ u p )+Vu p -($A p ) = 0, (2.4) 

where x p is the particle position, u p particle velocity, p p particle density, and V p particle 
volume. A p is the particle acceleration and F p = m p A p the total force acting on the 
particle of mass m p and are given below. Here, y X ' and Vu p ' are the divergence operators 
with respect to space and velocity, respectively. The individual particle positions and 
velocities can be obtained by solving the Liouville equation in Lagrangian framework for 
each particle p: 


d 

dt 


( x p) 


(2.5) 


p dt 


( u p) — Pp 


( 2 . 6 ) 


2.2.1. Particle Forces 

The main issue is to model the force on a particle. This may consist of the standard 
hydrodynamic drag force, dynamic pressure gradient, gradient of viscous stress in the fluid 
phase, a generalized buoyancy force, and inter-particle collision. The total acceleration 
of the particle A p is given as, 


A p — D p (u f — u p ) VPp+ ( 1 ) g + B p + A cp (2.7) 

Pp \ Pp/ 

Here B p is the generalized buoyancy force and A cp is the acceleration due to inter- 
particle forces. If Pp » p g the pressure gradient, viscous, and buoyancy terms are 
usually negligible. In the present study, the generalized buoyancy force is also neglected 
for simplicity. It is shown later that, even without the presence of this buoyancy force, 
one can obtain lift of particles in a shear flow. The drag force is caused by the motion of 
a particle through the gas. The standard expression for D p is used 


(2 - 8) 

O Pp JXp 

where Cd is the drag coefficient and is given by (Gidaspow 1994; Andrews & O’Rourke 
1996) 


C d = — (l + aReDQj 1 - 8 , for Re p < 1000 (2.9) 

= O.440J 18 , for Re p > 1000 (2.10) 

where C d is the drag coefficient for spherical particles, R p = {SVp/^Tt) 1 ^ is the particle 
radius. The particle Reynolds number ( Re p ) is given as 

Re p = 2 Pf & f\ u f~ u p\ R P . ( 2 . 11 ) 

df 

There is an indirect collective effect in this drag term: when there is a dense collection 
of particles passing through the fluid interphase momentum exchange term in equation 
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(2.2) will cause u ff to approach the particle velocity, u p thus decreasing the drag on a 
particle, a drafting effect. 

The probability function 4> integrated over velocity and mass gives the probable number 
of particles per unit volume at x and t in the interval (u p , u p + du p ) , (p p , p p + dp p ) , 
(V p , V p + dV p ). The particle volume fraction (0 p ) is defined from the particle distribution 
function as, 


0 P = 


QVpdVpdppdUp. 


(2.12) 


From continuity, the gas-phase volume fraction is obtained as 0/ = 1— 0 p . The interphase 
momentum transfer function per unit volume in equation (2.2) is given as 


F = 




D p (u f 


Up) 


1 

— VPp 

Pp 


dVpdppdUp 


(2.13) 


2.2.2. Collision Force 

The acceleration of particles due to inter-particle interactions (A cp ) is an important 
term in dense particulate and two-phase flows. For dilute and lightly loaded configura- 
tions, the particle volume fraction (0 P ) is small (< 10%), the inter-particle collisions are 
negligible and probability of particles overlapping each other is low. For dense partic- 
ulate flows, however, the particle volume fraction should not exceed the close-packing 
limit (which is usually around 0.6 for three-dimensional case). In the Eulerian-Eulerian 
approach for two-phase flows, this is ensured by force due to the gradient of interparti- 
cle stress in the averaged momentum equation for the particle phase (Gidaspow 1986; 
Gidaspow 1994). Same model was used in Eulerian-Lagrangian approach by Andrews & 
O’Rourke (1996), Patankar & Joseph (2001b), Snider et al. (1998), and Snider (2001). 
Accordingly, the expression for acceleration of particles due to collision (A cp ) is 


A — — 

Jr ^-cp — 


1 


OpPp 


VT 


(2.14) 


where r is the interparticle stress that provides a pressure type force that prevents packing 
of particles beyond the close-packing limit. Expression for r is given as (from Snider et 
al. (1998)) 


PsO 0 p 
O cp Op 


(2.15) 


where P s has units of pressure, 0 cp is the particle volume fraction at close packing and /? is 
a constant. The values for P s and /3 are obtained from Snider et al. (1998). In this model, 
the inter-particle acceleration due to particle collision is assumed to be independent 
of its size and velocity. This model is least expensive in terms of computational cost 
as particle binary pairs are not formed and the collision force is directly obtained by 
interpolation from equation (2.14). This model, however, couldn’t completely prevent the 
particle volume fraction to exceed the close packing limit and some numerical instabilities 
were experienced in the present unstructured code. Further research on this model will 
be conducted to eliminate this problem. An alternative but expensive collision scheme 
based on the distinct element method (DEM) of Cundall & Strack (1979). This scheme 
can be readily applied to parcel techniques used in Lagrangian spray simulations. A 
parcel represents a number (N p ) of droplets/particles of same size, velocity, and other 
properties. The effective radius ( R par ) of a parcel based on its total volume is then given 
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by 


R 


par — 


3 N P V P 

47T 


1/3 


(2.16) 


In this method, the interaction among particles and between wall and particles is taken 
into account separately to ensure that 0 P < 0 cp . The force F p ~ p on parcel p due to 
collision with parcel j is given by 


Fpj 11 /hr dpj ' ( Rpar ■ P T Rpan j T cf) 

(j^cdpj ejc (iiy, u, j • n pj'j n pj for d P j <C ( Rpar - P A Rpanj T o ) 

S p j (RpanP~k Rpar ■ j A Oi) d p j 

F p ~ p = -F p ~ p 
jp pj 

where d P j is the distance between the center of the p th and j th parcels, n p j is the unit 
vector from the center of parcel j to that of parcel p, a is the force range, k c the stiff- 
ness parameter, and rj c the damping parameter. Tsuji et al. (1993) used the following 
expressions to compute the damping parameter 


Vc = 2 a 


a = — 


/ m p k c 
In (e/n) 


where e is the coefficient of restitution. Similarly, the parcel-wall force ( Fi p w w ) on parcel 
p due to collision with wall w is 


f p~ w = 

pw 


fipw — 


0 for dp W Z/ ( RpariP ^) 
(k c fij w 7/ c (iip) • n /r „, n pw 
(Rpar i P T oO d pw 


for 


dpw ^ ( RpariP (f) 


where d pw is the distance between the parcel and the wall, and n pw is the unit vector from 
the wall to the center of the parcel. The total collision force is obtained by looping over 
all particles and walls. The corresponding particle acceleration is obtained by dividing 
the collision force by parcel mass (m p = Npp p V p ). 

3. Numerical Method 

The collocated, finite-volume numerical scheme on unstructured grids developed by 
Mahesh et al. (2003) is modified to take into account the gas-phase volume fraction. In 
addition, the particle centroids are tracked using the Lagrangian framework developed by 
Apte et al. (2003b). The important feature of the numerical scheme is the computation 
of 0/ and 0 P on the unstructured grids. As the particles in the simulations performed 
do not move out of the domain or are destroyed, the total volume occupied by the 
particles remains constant from mass-conservation. As a particle crosses a particular 
control volume, its contribution to the particle volume fraction on neighboring cells must 
provide global mass-conservation. In addition, it was observed that the Eulerian volume 
fraction field should also be smooth in order to avoid numerical instabilities encountered 
due to changes in 0 as the particles move from one grid cell to another. A strategy similar 



166 


Apte, Mahesh, & Lundgren 


Computational domain, 0.2 x 0.6 x 0.0275m 
Fluid density, 1.254 kg/m 3 
Numer of Parcels, 2880 
Diameter of particles, 500 gm 


Grid, 10 x 30 x 5 
Particle Density, 2500 kg/m 3 
Particles per parcel, 3375 
Initial particle concentration, 0.4 


Table 1 . Parameter description for the gravity-dominated sedimentation case. 



Figure 2. Temporal evolution particle distribution during gravity-dominated sedimentation. 
Initially the particles are randomly distributed over the box rectangular box. 

to the two-way coupling methodology by Maxey & Patel (2001) is used to compute the 
volume fraction. The interphase momentum transport terms are also treated in a similar 
way. 

The particle equations are integrated using third-order Runge-Kutta schemes for ode- 
solvers. At each Runge-Kutta step, the particles were located and the collision force 
was re-computed. This was found necessary to ensure that the close-packing limit for 
particle-volume fraction is not exceeeded. 


4. Results 

We first simulated sedimentation of solid particles under gravity in a rectangular box. 
Details of this case are given in Table 1. The initial parcel positions are generated ran- 
domly in the top half of the box, 0.3 < y < 0.6. These particles are then allowed to settle 
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Computational domain, 0.2 x 0.6 x 0.0275m 

Gas jet velocity, 9 m/s 

Fluid density, 1.254fcg/m 3 

Numer of Parcels, 2880 

Diameter of particles, 500 gm 


Grid, 10 x 30 x 5 
Jet diameter, 0.04m 
Particle Density, 2500kg /m 3 
Particles per parcel, 3375 
Initial particle concentration, 0.4 


Table 2. Parameter description for the simulation of fluidization by a gas jet. 


through the gas-medium under gravity. The collision model is important here near the 
bottom wall after the initial group of parcels hit and bounce back from the bottom wall. 
This prevents particle volume fraction to exceed the close-pack limit. The upper mixture 
interface between the particles and the fluid is closely approximated by h = gt 2 / 2 similar 
to that obtained by Snider (2001). The particles eventually settle down with close-packing 
near the bottom wall. Snider (2001) and Patankar & Joseph (2001a) did sedimentation 
under gravity, in an inclined pipe with large number of parcels and used the particle 
stress model for collision force. This allowed them to use large number of parcels as the 
computational cost for this collision model is negligible. We plan to use this model to do 
similar validation test cases using the unstructed LES code. One difficulty at the time 
of writing was that this model didn’t always prevent the particle volume fraction from 
exceeding the close-pack limit. Further investigations on this matter will be performed. 

4.1. Gas-Solid Fluidization 

We consider the problem of fluidization of solid particles arranged in an array at the 
bottom of a rectangular box. Fluidization is achieved by a jet of gas issuing from the 
bottom of the box. The flow parameters are given in Table 2. The choice of collision 
parameters is based on the recent computations by Patankar & Joseph (2001b). The 
particle motion is mostly dominated by the hydrodynamic drag force and collision model 
should not affect the overall particle motion. The collision model, however, is important 
in governing the particle behavior near the walls and helps prevent the volume fraction 
from exceeding the close-pack limit. 

Figure 3 shows the position of parcels at different times during bubbling fluidization. 
Parcel diameters are drawn to scale. The jet issuing from the bottom wall pushes the 
particles away from the center region and creates a gas-bubble in the center. The particles 
collide with each other and the wall and are pushed back towards the central jet along the 
bottom wall. They are then entrained by the jet and are levitated. This eventually divides 
the central bubble and two bubbles are trapped. The particles tend to move upwards and 
collide with the upper wall and remain levitated during future times. The computational 
results are in good agreement with the simulations of Patankar & Joseph (2001b) as 
well as in qualitative agreement with experiments on jet fluidization. Similar results 
are reported using Eulerian-Eulerian approach in two-dimensions by Ding & Gidaspow 
(1990). 


4.2. Fluidization by Lift of Spherical Particles 

The transport of particles by fluids in coal-water slurries, hydraulically fractured rocks 
in oil-bearing reservoirs, bedload transport in rivers and canals and their overall effect 
on the river bed erosion etc., are important scientific and industrial issues in particulate 
flows. In order to understand fluidization/sedimentation in such conduits, (Choi & Joseph 
2001) performed a DNS study of fluidization of circular cylinders (300 particles) arranged 
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Figure 3. Temporal evolution of particle distribution during fluidization by a gas jet. Initially 
all the particles are uniformly arranged in layers at the bottom of the rectangular box. Air 
is injected through a rectangular slot at the bottom wall. Air bubbles are trapped within the 
particles and the growth and pattern of these bubbles are in agreement with experimental 
observations. 


Computational domain, 63 x 12 x 12cm 

Gas jet velocity, 9 m/s 

Fluid density, 1 g/cm 3 

Particle Density, 1.01g/cm 3 

Numer of Parcels, 3780 

Initial array height, 4.75cm, 

Pressure gradient, 20 dyne /cm 3 


Grid, 20 x 11 x 10 

Jet diameter, 0.04m 

Fluid viscosity, 1 poise 

Diameter of particles, 0.95cm 

Particles per parcel, 1 

Initial centerline velocity, 360 cm/s 


Table 3. Parameter description for the simulation of fluidization of spherical particles in a 

plane Poiseuille flow. 


at the bottom of a channel in plane Poiseuille flow. They observed that with sufficient 
pressure gradient across the channel, the particles initially at rest in the lower half of the 
channel start moving and roll over the wall. Particle rotation in a shear flow generates lift 
and the channel is fluidized after some time. We attempt to capture this phenomenon by 
our Eulerian-Lagrangian model. The flow parameters are given in Table 3. As opposed to 
Choi & Joseph (2001), we are performing three-dimensional simulations and our particles 
are spheres. As shown in Fig. 4 we observe that the effect of volume displacement due 
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Figure 4. Temporal evolution of particle distribution during fluidization by lift in a plane 
Poiseuille flow. Also shown are contours of axial gas-phase velocity. Initially all particles are 
arranged at the bottom of the channel. 


to particles is to set up Kelvin-Helmholtz type waves between the fluid and particle 
layers. A vertical pressure gradient is created and gives vertical velocity to the particles 
and the channel gets fluidized. We also did several test cases, with higher grid resolution, 
increased density ratios to obtain similar results. With increased particle density, the two- 
way coupling between the particle and fluid momentum equation, decelerates the fluid in 
the bottom half of a channel and an inflection point is created in the axial velocity profile. 
This eventually cause lift and particle dispersal. In liquid-fuel combustor applications, 
we believe that similar mechanism can be observed in the dense-spray regime near the 
injector and will allow us to capture the important features of primary atomization. 


5. Summary 

We have developed a numerical scheme which accounts for the displacement of the fluid 
by particles and is applicable to dense particulate flows and spray regimes close to the 
nozzle injector. This scheme has been implemented into the unstructured LES code. We 
performed various numerical simulations of fluidization and particulate flows to validate 
our scheme. It was shown through simulation of channel flow with layers of spherical 
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particles at the bottom wall that the effect of particle volume fraction is to displace the 
fluid in a control volume and create vertical pressure gradient through two-way coupling 
between the particles and fluid. This gives Kelvin-Helmholtz type instability waves which 
result in dispersal and lift-off of particles. This liftoff of dense-particles cannot be observed 
without explicit modeling of lift-forces on particles by using the standard point-particle 
approximation in Lagrangian-Eulerian framework. Similar K-H waves are present in the 
dense spray regime near the injector nozzle which may affect the overall spray dynamics 
and breakup in practical combustors. We plan to apply this model coupled with secondary 
breakup model to perform simulations of spray atomization. We also plan to perform 
several validation studies of this methodology for dense particulate/granular flows to 
address some issues related to collision scheme, efficient and consistent computation of 
particle volume fraction on multiple domain-decomposition. 
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